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Abstract 

Rapid retreat of ice in the Amundsen Sea sector of West Antarctica may cause drastic sea 
level rise, posing significant risks to populations in low-lying coastal regions. Calibration 
of computer models representing the behavior of the West Antarctic Ice Sheet is key for 
informative projections of future sea level rise. However, both the relevant observations 
and the model output are high-dimensional binary spatial data; existing computer model 
calibration methods are unable to handle such data. Here we present a novel calibration 
method for computer models whose output is in the form of binary spatial data. To mitigate 
the computational and inferential challenges posed by our approach, we apply a generalized 
principal component based dimension reduction method. To demonstrate the utility of our 
method, we calibrate the PSU3D-ICE model by comparing the output from a 499-nrember 
perturbed-paranreter ensemble with observations from the Amundsen Sea sector of the ice 
sheet. Our methods help rigorously characterize the parameter uncertainty even in the 
presence of systematic data-model discrepancies and dependence in the errors. Our method 
also helps inform environmental risk analyses by contributing to improved projections of sea 
level rise from the ice sheets. 

Keywords: Computer Experiments, Spatial generalized linear mixed models, Climate change, 
Gaussian processes, Principal components 

1 Introduction 

Mass loss from the polar ice sheets has the potential to make substantial contributions to future 
sea level rise. The ice sheets therefore pose substantial risks to people living near present sea 
level. Taken together, the ice sheets contain sufficient ice to raise global mean sea level by 
up to 60 meters if they were to melt completely; the West Antarctic Ice Sheet (WAIS) might 
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contribute up to 4 meters to future sea level rise (Fretwell et al. 2013), whereas the Greenland 
and East Antarctic Ice Sheets might contribute up to 7 meters and up to 50 meters to sea level 
rise, respectively. Recent observations of increased surface air temperatures at high northern 
latitudes, and warmer ocean waters at high southern latitudes, suggest that some fraction of 


this ice may melt over the next few decades or centuries (Zhang 2007; Serreze and Barry 2011). 
Geologic observations also indicate that the ice sheets can make rapid contributions to sea level 


(e.g., Deschamps et al. 2012). A significant fraction of the world’s population lives near present 


sea level (Nicholls et al., 2008), and these people will face an increased risk from storm surges 
as sea level rise progresses. The future behavior of the ice sheets is hence a key input for risk 


analyses (e.g., Lempert et al. 2012). 

Despite the ice sheets’ importance for future sea level rise, projecting their future behavior 
using computer models remains challenging. In contrast to global climate models, which repre¬ 
sent the behavior of the atmosphere, oceans, vegetation, and other parts of the Earth system, 
ice sheet models typically include advanced treatments of ice flow and simplified representations 
of the boundary conditions at the ice sheets’ interfaces with the atmosphere and oceans. Until 
recently, many ice sheet models lacked realistic representations of processes that may play an 


important role in the ice sheets’ future development (e.g., Little et al. 2007). Recent studies 


(see, e.g. Stone et al. 2010 Applegate et al. 2012) confirm that imperfect knowledge of model 
input parameters leads to large uncertainties in modeled sea level contributions. Thus, char¬ 
acterizing and reducing parametric uncertainty is a key step in producing ice sheet projections 
that can properly inform climate risk analyses. 

Calibration of ice sheet models is challenging because the relevant data are high-dimensional 
binary spatial data, and the models themselves are computationally expensive. Recent efforts 


to calibrate ice sheet models include Gladstone et al. (2012), McNeall et al. (2013), and Chang 


et al. (2014a). These studies represent significant advances in the state of the art for ice sheet 


model calibration, but each has significant limitations. Gladstone et al. (2012) uses a heuristic 
approach which chooses input parameter settings that yield output values which fall in the 


95% confidence set. The confidence set is defined by the ‘three-sigma rule’ (Pukelsheim, 1994), 
choosing the margin of errors for the confidence sets as ±3 standard errors, and a Bonferroni- 
type multiple test procedure to generalize it to the multivariate case. The approach ensures 
robustness against assumptions about the model-observation discrepancy and the observational 
errors. The framework requires running the ice sheet model for at least thousands of parameter 
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settings and hence is applicable only to heavily simplified ice sheet models such as the one¬ 


dimensional flow-line model used by Gladstone et al. (2012). McNeall et al. (2013) uses a 
Gaussian process emulator for three quantities: the overall ice volume, the maximum ice height, 
and the area covered by ice. The use of Gaussian process emulators makes this approach 
computationally feasible, even with computationally expensive 3-D ice-sheet models. However, 
their analysis is based on heavily aggregated quantities and therefore does not utilize spatially 
resolved information. Moreover, instead of utilizing the probability model given by the emulator, 
they rely on a heuristic approach based on an “implausibility measure”, which is the model 


output mean standardized using the observational data. The approach in Chang et al. (2014a) 
improves upon these two approaches by using the probabilistic calibration approach proposed by 


Kennedy and O’Hagan (2001) and extended and generalized by many others (see, e.g. Kennedy 


and O’Hagan 

2001 

Bayarri et al. 

2007 

Higdon et al. 

2008) 


However, current calibration approaches are unable to utilize important sources of infor¬ 
mation about ice sheet behavior such as the extent of ice covered area in modern- and paleo- 


observations (Alley et al., 2010) due to the binary nature of the model output and observations. 
For example, some Greenland ice cores contain ice from the last interglacial period (from 125,000 
years ago), whereas others do not; however, rigorous calibration approaches that can take this 


type of data into account are lacking. Chang et al. (2014a) circumvents this problem by us¬ 
ing datasets obtained by aggregating the 2-dimensional spatial pattern of ice thickness into a 
1-dimensional profile to make the data more suitable for Gaussian process-based approaches 
and hence is silent on using datasets that are in the form of non-Gaussian spatial processes. 
Aggregating data in this manner likely results in significantly larger projection uncertainties 


than would result if unaggregated data were used (Chang et al., 2014b). Thus, more work to 
properly incorporate binary spatial data into ice sheet model calibration is needed. 

Here we propose a new calibration approach for high-dimensional binary spatial data, ap¬ 
propriate for calibration of ice sheet models using data on the spatial pattern of ice coverage. 
The spatial data that we use here are gridded binary patterns, where the outcome at each lo¬ 
cation represents absence or presence of grounded ice. Our novel calibration approach extends 
the Gaussian process-based calibration framework to high-dimensional binary spatial data. Our 
approach is based on a generalized linear model framework with a latent Gaussian process. To 
mitigate the considerable computational challenges posed by the high-dimensional latent pro¬ 
cess, we reduce the dimensionality of the model output and the data using a logistic principal 
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component analysis (Collins et al. 2002 Lee et al. 2010). This approach avoids inference with 
high-dimensional latent variables while also sidestepping the computational burden posed by 
high-dimensional matrix computations. 

We use this approach to calibrate an innovative ice sheet model (PSU3D-ICE; see below), 
using data on the spatial extent of grounded ice (i.e., binary pattern of presence and absence of 
grounded ice) within the Amundsen Sea sector of the West Antarctic Ice Sheet. The Amund¬ 
sen Sea sector is likely the largest contributor to sea level rise of any Antarctic ice drainage 
(Pritchard et ahj 2012). It includes the rapidly thinning Pine Island and Thwaites Glaciers, 
which drain ice from the interior of the West Antarctic Ice Sheet into the ocean. The glaciers 
terminate in ice shelves where the grounded ice becomes too thin to remain in contact with 
the seafloor. Ice velocities are high within the glacier trunks, but slow elsewhere. The position 
of the grounding line, or the narrow region separating the grounded ice from the floating ice 


shelves, is known from modern geophysical measurements (Fretwell et al., 2013). 

The Pennsylvania State University three-dimensional ice sheet model, PSU3D-ICE, simu¬ 
lates the combined ice sheet-ice shelf system using a so-called “hybrid” dynamical core, which 
smoothly merges the shallow-shelf and shallow-ice approximations to the Stokes equations de¬ 


scribing ice flow (Kirchner et al., 2011). This approach allows the model to simulate the full 
ice sheet-ice shelf system realistically, while minimizing the computational cost of the model. 
Pollard and DeConto (2009, 2012bfa) provide detailed descriptions of the model, which has been 


part of several ice-sheet model intercomparisons (e.g., seaRISE, Bindschadler et al., 2013). 

The remainder of this paper is organized as follows. In Section [2] we describe the model 
output and the observational data that we use in our analysis. In Section [3] we introduce the 
basic framework for calibration using binary spatial data and explain the computational chal¬ 
lenges posed by high-dimensional spatial data. We formulate a reduced-dimensional approach 
to mitigate these challenges in Section [4} We describe how to compute the logistic principal 
components (PCs) for binary computer model output and formulate an emulation framework 
based on them, and then set up a calibration model that relies on the reduced-dimensional 
emulator. Section [5] describes the implementation details and results for our experiments on a 
simulated example as well as the application of our methods to the ice sheet model for WAIS 
described in Section [2] We also discuss the preliminary implications of our research on the 
study of the behavior of the West Antarctic Ice Sheet. We conclude this paper with a summary 
and future directions in Section [6l 
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2 Model Description and Observational Data 


We use an innovative 3-D ice-sheet model, which uses boundary-layer techniques to capture chal¬ 
lenging grounding-zone dynamics while still being efficient enough to feasibly perform long-term 
simulations. Basically the model simulates the long-term evolution of ice extent and thickness, 
changing due to slow deforming flow and basal sliding under its own weight, surface snowfall and 
melt, ocean melting under floating ice shelves, and iceberg calving, over the Antarctic continent. 
For the simulations in this manuscript, the model is applied to a nested domain spanning West 
Antarctica, with a relatively fine grid resolution of 20 km. This nested grid is driven at its 
lateral boundaries by stored results from a previous coarser-scale continental simulation. 

The model is run over the last 40,000 years, initialized appropriately at 40 ka (40,000 years 
Before Present or BP, relative to 1950 AD) from a previous long-term run. Atmospheric forcing 


is supplied using a modern climatological Antarctic dataset (ALBMAP: Le Brocq, A.M. and 


Payne, A.J. and Vieli, A)| 2010), with uniform cooling perturbations applied proportional to 


a deep-sea-core <5 ls O record (Pollard and DeConto, 2009, 2012a). Oceanic forcing is supplied 


using archived ocean temperatures from a coupled AOGCM simulation of the last 20 kyr (1 


kyr=l,000 years) (Liu et ah, 2009b). Sea level versus time, rising during this period primarily 


due to Northern Hemispheric ice-sheet melt, is prescribed from the ICE-5G dataset (Peltier 


2004). Modern bedrock elevations are obtained from the Bedmap2 dataset (Fretwell et al. 


2013). 


After reaching present day (0 ka), each run is extended for 5,000 years with a very basic rep¬ 
resentation of future warming. Atmospheric and oceanic temperatures are uniformly increased 
by 6 and 2 °C, respectively, ramped up linearly from 0 to 150 years after present and held 
constant thereafter. Ocean-temperature increases are confined to a longitudinal sector around 
the Amundsen Sea Embayment of West Antarctica, corresponding to the location of observed 
sub-ice melt increases in recent decades. This simple prescription of future temperature changes 
causes drastic ice retreat into the West Antarctic interior in many of the runs. More realistic 
future warming scenarios are planned for future work. Other aspects of the model are described 
Pollard and DeConto (2012bja). 


m 


All 3-D ice sheet models have a number of poorly constrained internal parameters. To 
evaluate the effects of these uncertain parameters in the PSU3D-ICE model, we selected 10 
highly uncertain parameters to vary in a Latin hypercube (LHC) ensemble (Table [I]). LHC 
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design (see, e.g. McKay et al. 1979 Stein, 1987) allows us to explore a reasonably wide range of 


parameter values with a relatively smaller number of model runs even with a high-dimensional 


parameter space. We use a random LHC design generated by ‘lhc’ package in R (Carnell, 2012) 


While a random LHC can result in a design with highly correlated input parameter values (cf. 


Joseph and Hung, 2008), we have confirmed that in our design the correlation coefficients 


between the parameter values are very low (all less than 0.08). Four of these parameters 
are generally recognized as particularly important, with strong effects on past and present ice 


configurations, but which represent uncertain processes in glaciology. See Section |5.2| for a 
description of these parameters. Each parameter value is remapped to the [0,1] interval with 


0.5 generally being the non-statistically tuned value by previous research (Pollard and DeConto 


2012a Pollard et al. 2015) and 499 design points for the input parameters are chosen by LHC 


design for the 10-dimensional unit cube. 

We calibrate the model by inferring these ten parameters using the observed modern day 
grounding-line geometry, i.e., the modern “coastline” map of grounded ice versus floating ice, 


from a modern Antarctic dataset (Bedmap2; Fretwell et al., 2013). The original model output 
covers the entire WAIS, but we use a subset of the output that corresponds to the observational 
data for ASE (Amundsen Sea Embayment). The spatial pattern for each parameter setting 
and the observational data has 86 x 37 regular grid points. At each grid cell, we have a binary 
outcome with 1 representing presence of grounded ice and 0 representing absence of grounded 
ice. See Figure [l] for the observational data and an example model output. Note that both 
the model output and the observational data are in the form of high-dimensional and binary 
spatial patterns; the challenges posed by these data structures will drive much of our methods 
development in the following sections. 


3 Computer Model Emulation and Calibration Using Binary 
Spatial Data 

In this section, we introduce a framework for computer model emulation and calibration with 
Gaussian spatial data and then describe our new approach for emulation-calibration with binary 


spatial data. We follow the two-stage approach to computer model calibration (Bayarri et al 


2007, Bhat et al., 2012 Chang et al. 2014b, 2015). In the two stage approach, the computer 


model is first approximated using a Gaussian process emulator, and then in the calibration step, 
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the model parameters are inferred based on observational data using the emulator. Separating 
the emulation step from the calibration step makes it easier to diagnose the performance of 
the emulator and improves the identihability of the parameters in the calibration model (see 


Rougier, 2008; Bhat et ah, 2010 2012, for more discussion). Here the computer model output 


is in the form of binary spatial data, and we specify a generalized linear model framework and 
formulate an emulator that interpolates natural parameters across the model input settings using 
a Gaussian process. For the calibration step we introduce a model for data-model discrepancy 
and model the natural parameters for the observational data using the constructed emulator 
and the discrepancy model. When the model output is in the form of Gaussian spatial data, 
there is a possibility that we underestimate the parametric uncertainty when we do not account 


for the uncertainty in the estimated emulator parameters (see for instance, Bayarri et al. 2007 


Bhat et al. 2012; Liu et al., 2009a; Chang et ah, 2014b). However, when the model output is 


in the form of binary spatial data, based on several simulated examples we find that we do not 
underestimate the parametric uncertainty. In fact, allowing uncertainties for this parameters 
leads to identihability issues in the inference in the calibration stage. 

Computation presents considerable challenges. We therefore propose a principal components- 
based framework in order to make our approach computationally feasible. To help set the stage 
for the description of our framework for emulation and calibration in this manuscript, we begin 


with a description of emulation-calibration using Gaussian spatial data in Section 3.1 This 


framework is also described in Kennedy and O’Hagan (2001); Bayarri et al. (2007); Higdon 


et al. (2008); Bhat et al. (2012); Chang et al. (2014b). Sections 3.2 and 3.3 provide details 


about our framework for emulation and calibration for binary spatial model output and obser¬ 
vational data. Section 3.4 describes the computational and inferential challenges, and Section [4] 
provides details about our computationally expedient reduced-dimension approach to mitigate 
these challenges. 

We use the following notation henceforth. Let H(0,s) be the computer model output at 
a parameter setting 6 G 0 and a spatial location s G S where 0 C R d is the parameter 
space and S is the spatial held that covers the geographical area of interest. The integer d 
represents the number of input parameters that the computer model has. We assume that we 
have obtained model runs at p design points 0\,... ,0 p G 0 and each model run has n spatial 
locations si, ..., s n G S. We let Y be a p x n matrix of the computer model output where the 
(*, j)th element Y(6i, s j) is the model output at the ith parameter setting 6i and the jth spatial 
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location s j. Similarly, let Z be the n-dimensional observational data where its jth element Z(sj) 
is the observation at the j th spatial location s j. For the ice sheet model calibration example in 


Section 5.2, p = 499 and n = 86 x 37 = 3,182. 


3.1 Computer Model Emulation and Calibration Using Gaussian Spatial 
Data 


We briefly outline emulation and calibration for Gaussian spatial data (cf. Bhat et ah, 2010 


2012, Chang et al., 2014b 2015). The purpose of this subsection is to make it easier for readers 


to follow the remainder of this section. Here we assume that Y(8, s) can be reasonably modeled 
by a Gaussian distribution. In the emulation stage, we fit the following Gaussian process model 

^ ^ >-0r. 

and the spatial locations si,..., s n simultaneously: 


for the computer model output Y to interpolate the values at the parameter settings 0 \,..., w p 


vec(Y)~X(X/3,S(£ y )), (1) 

where vec(-) is the vectorization operator that concatenates columns into one vector. The rip x b 
matrix X contains values of b different covariates that can be used to represent general trends 
in the model output. The np x np covariance matrix X(£ y ) with a vector of parameters is 


and spatial locations si,..., s n . In general, the mean term X/3 can be set to 0 with a proper 
centering of Y, and the covariance parameters in £ y can be estimated by the maximum likelihood 
estimate £ y . 

In the calibration step, we model the observational data using the following model: 


defined by a positive definite covariance function for a Gaussian process (see, e.g. Sacks et al. 


1989) that interpolates between the model output at different parameter settings 


Z = t7(0,Y) + <5, 

where r](0, Y) is the n-dimensional emulated vector, a computer model output approximated 
by interpolating the model runs in Y using the model in ([!]), and <5 r\_/ N(0, £(£ d )) is an n- 
dimensional spatially correlated zero-mean Gaussian vector with covariance parameter that 
represents the systematic model-observation discrepancy. By choosing standard priors for the 
parameters in the model, one can infer 6 via Markov Chain Monte Carlo (MCMC) along with 
other parameters. 
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3.2 Computer Model Emulation Using Binary Spatial Model Output 

In this section and for the remainder of this paper, we assume that Y(9, s) is binary with 0-1 
values. We assume the model output Y(9i,Sj ) to be a Bernoulli random variable with success 
probability p % j = P(Y(9i,Sj ) = 1). We denote the corresponding logit log ^^ by 7 ij and 
hence the probability mass function for Y(9j,Sj ) can be written as 


P{Y(0 i ,s j ) = y ij )=p™(l-p ij ) 1 -™ 

= ( ex P ( 70 -) v Mi ( 1 

V 1 + exp ( 7 ^)) \ 1 + exp( 7 ij )) 

=(1 + exp(-(2yij - 1 ) 7 ij)) -1 . 

We let r be the pxn matrix of natural parameters for model output, where its (i, j)th element is 
7 ij. By assuming that the elements in the model output Y are conditionally independent given 
the natural parameters in T, the likelihood function for the model output Y can be written as 
follows; 


p n 


p n 


L(Y|r) = n n P(Y(9i, Sj) = y l3 ) = + exp(-(2 yii - l)^))" 1 . (2) 

i=l j=l i= l j=l 

The dependence among the elements in Y is modeled through the dependence among the 
elements in T. To approximate the computer model we interpolate the natural parameters at 


different input parameter settings using a Gaussian process (cf. Sacks et ah, 1989) with zero 


mean and a covariance function that depends on the input parameter settings and the spatial 
locations for model output, i.e., 


Cov { n /ij,Jkl) = C (9i,9 k , Sj , sr ,£ y ) , 


(3) 


where C is a valid covariance function over (0x5 ) 2 and is a vector of covariance parameters. 
Note that estimating the natural parameters {7 ij,i = 1... ,p, j = 1,... ,n} is not possible by 
simply maximizing Q since this is an ill-posed problem without any constraints on them (see 


Section 3.4 for details). If estimates of the natural parameters are available, we can fit the Gaus¬ 
sian process to those natural parameters by maximizing the likelihood function corresponding 
to the probability model 

vec(r)~JV(0,E(£ y )) 
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with respect to £ , where the covariance matrix £(£) is given by a function of the form in ([3]). 
We then use the fitted process to predict the vector of natural parameters at any new parameter 
value 9. We call the fitted Gaussian process an “emulator” and the resulting predictive process 
an “emulator process”. 

3.3 Computer Model Calibration Using Binary Spatial Data 

Again, in this section and for the remainder of this paper, we assume that Z( s) is binary with 
0-1 values. We set up a calibration model that links the input parameter to the observational 
data using the emulator constructed in the previous section while considering the systematic 
discrepancy. Let A = (Ai,..., \ n ) T be the n-dimensional vector of natural parameters for the 
observational data, then our calibration model is given by 

A = r?(0,Y) + <5, (4) 

where 8 is an n-dimensional vector for the model-observation discrepancy. Here we can think 
of 9 as a “fitted” value for the parameter 9, that is, 9 at which the values of the natural 
parameters match those from which the observations are assumed to originate. For 8 we use 
an n-dimensional Gaussian random vector with covariance parameters that is independent 
of the emulator process r/ (9, Y). This corresponds to the following probability model for each 
observation given by the standard logistic regression framework: 

= / exp(Aj) y ( 1 A 1 ~ Zi 

V 1 + exp(Aj) ) \ 1 + exp(Aj) / 

=(1 + exp(—(2 zj - l)A i )) _1 , 

where A j is the jth element of A. The log-likelihood function for A, 9 , and given T and Z 
can be written as 



e(z,r\\,e,z d ) = e(r\\,o,s d ) + e(z\\). 

Here ^(T|A, 9, £ d ) and £(Z|A) are given by the probability models in Q and ([5]) respectively. 
With a standard prior specification for the input parameters 9, the discrepancy parameters in 
and the latent process A, we can define the posterior density for those parameters and carry 
out Bayesian inference using MCMC. 
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3.4 Computational Challenges 


The framework described above faces computational and inferential challenges when the number 
of spatial locations n is large. The emulation step described in Section [3.2| requires estimation 
of np natural parameters. This is an ill-posed problem since the number of data points in the 
model output is also np. It is not straightforward to find a set of constraints that leads to well- 
posed parameter inference while still allowing enough flexibility for the resulting process. Even 
if estimates of the natural parameters are available, emulation still poses significant compu¬ 
tational challenges. Emulation involves numerically finding the maximum likelihood estimate 
of £ y , requiring repeated evaluation of a computationally infeasible likelihood function; each 
likelihood evaluation here involves Cholesky decomposition of the np x np covariance matrix 
whose computation scales as 0(n 3 p 3 ). For the example in Section [ 5 J this calculation trans¬ 
lates to a computational cost of | x p 3 x n 3 = | x 499 3 X 3,182 3 = 1.33 x 10 18 flops, which 
in turn corresponds to at least hundreds of thousands of hours of computing time for a 3.0 
Ghz single core. Moreover, storing np X np double-precision covariance matrix will require 
8 n 2 p 2 = 8 x 499 2 X 3,182 2 = 20,169.33 terabytes of memory space. The calibration step is sub¬ 
ject to similar challenges. Estimating the input parameters 6 and the discrepancy parameters 
requires integrating out the n-dimensional discrepancy process 5, which makes it difficult to con¬ 
struct a well-mixed MCMC algorithm even for a moderate size of n. Moreover, each likelihood 
evaluation involves Cholesky decomposition of the n x n covariance matrix and its computa¬ 
tional cost scales as 0(n 3 ). This translates to a computational cost of ^ x 3,182 3 = 3.59 x 10 10 
flops in the example in Section [5j which corresponds to 30 minutes of computing time for a 3.0 
Ghz single core. 

4 Reduced-Dimension Approach for Binary Spatial Data 

We propose a reduced-dimension approach to overcome the computational challenges discussed 
in the previous section. In the emulation step, we build an emulator for the logistic principal 
components of the model output to solve the inferential issues in estimating np different natural 
parameters and the computational issues for dealing with the np x np covariance matrix. Emu¬ 
lation based on principal components helps us to reduce information loss caused by dimension 
reduction by focusing on the principal directions for the natural parameters that have the max¬ 
imum variation across the input parameter settings 6 1 ,..., 6 P . In the calibration step, we use 
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a basis representation approach to model the discrepancy process, which removes the need to 
integrate out the n-dimensional discrepancy process and compute the Cholesky decomposition 
of the n x n covariance matrix. 


4.1 Dimension Reduction for Binary Model Output 


In this subsection, we describe the logistic principal component analysis for binary model output 


in detail by closely following Lee et al. (2010). For dimension reduction we assume that T can 
be written as 


r = lpn 1 


WK 


T 

y > 


( 6 ) 


where K y is an n x J y orthogonal basis matrix, W is the p x J y principal component matrix 
with (i,j )th element ), and /x is the n X 1 mean vector. We rewrite the likelihood in ([2]) 
as 

p n 

~ (7) 


L(Y|p, K y , W) = []n 9(y*iM + w * k L'))’ 

i=lj=l 

where g(x) = (1 + exp(—x))^ 1 , pj is the jth element of p, y*- = 2y^ — 1 , w, is the ith row of 


W, and k ;(y j is the jth row of K 


y 


Following Lee et al. (2010), we estimate the matrix W and K„ by maximizing the likelihood 


in 0 using the majorization-minimization (MM) algorithm ( |Lange et al.[ |2000| | Hunter and| 
L ange[ [2004). For each iteration of the MM algorithm we minimize the majorizing function 
instead of the negative likelihood function. The majorizing function is a function that is (i) 
of the same value as the negative likelihood at the current point of the iteration, (ii) greater 
than the negative likelihood for all other values of the parameters, and (iii) in an easier form 
to minimize than the original negative likelihood, such as a quadratic function. As a first step, 
we note that 


-logy(.T) < - log g{y) - (1 - g(y))(x - y) + -(x - yf 


( 8 ) 


for any value of y. The equality holds if x = y. The proof of this inequality can be found in 


de Leeuw (2006) (Theorem 4). We can complete the square by adding a term that does not 
depend on x to rewrite the right hand side as follows; 


1 2 

-log g{x) < - log g(y) + - (x - y -4(1 - g{y))) . 


(9) 


By the fact that can take only the value of either +1 or -1, plugging in y\fiij and y\j r i\- 


M 
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into x and y respectively yields 


-log9(^'7u) < - + ^ ( 


1 / (in) 

7ij ~ X ij 


( 10 ) 


where 7 ij is the ( i, j)th value of T, 7 ^" is the value of 7 y at the mth iteration, and = 7 + 
4y*-(l -giVij 7i, m) ))- Note that Q implies that 7 ^ = pj + w;k^- and 7 ^ m) = /x( m) + w,) m) k^ T 
where k|"j\ and w- m) are the values of pj, k yj , and w* at the mth iteration respec¬ 

tively. Since — log g(y*j7ij) does not depend on 7 ^, the majoring function for the negative log 
likelihood function in ([T]) is given by 


p n 


2 (4 m - /u - Wik L 

i= 1 j=l 


For each iteration the minimum of the majorization function occurs at 


Pj — 


1 p 

1 \ " ( (m) , T 

- I - Wjk, 


P 
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y,3 


i =1 

T1 


w i =(K y K y ) _i K y (x.\ — At), 


_ 1 K y (x 4 - m)T 
k,7 =(W T WpW T (x^ - l p /; y ) , 


( 11 ) 

( 12 ) 

(13) 


where and xP are the ztli row and the jth column of X( m ) respectively and X^ m ) is a 
p x n matrix having x^ as its (i, j)th element. Cycling through 0 - pi ) results in a local 
maxima of the likelihood function in Q. Results related to convergence for this algorithm can 
be found in Hunter and Lange (2004). Of course, as in other optimization algorithms one needs 
to be careful about possible stationary points. We summarize the algorithm as follows. 

1. Initialization: Set m = 1 and choose starting points ,..., fin j T , = 

(w[ 1)T , • • •, w« T ) and K y 1} = (k^J T , ..., k$T) • 

2. Compute 7 ^ m) = juj m) + wJ m) k^ )T , x = 7 ^ m) + 4y*(l - vppf” 0 )) and let X<» = 

K m) }- 

3. Update p m+1 ) = 1 (xM - W^K^Y l p . 


4. Compute W* = (X( m ) - (k^K^ 

decomposition W* = UAV T . Update W( m+1 ) = UV T . 

5. Update K { ™ +1) = (X( m ) - l p ^ m+1 '> T ) T W^ m+1 \ 


-1 


and find its singular value 
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6. Set m + 1 and repeat from 2 until convergence. 


4.2 Emulation Using Logistic PCs 

We fit a Gaussian process for each column of W separately to construct a ./^-dimensional 
emulator process. We denote the value of the jth principal component as Wj(0) and set up a 
Gaussian process for each jth column with the following exponential covariance function for 
any parameter settings 0 and O': 


Cov ( wj{6 ), Wj(O')) = Kj exp — 




2—1 


+ (, 1(0 = 0 ’), 


(14) 


"13 


with the partial sill Kj, the range parameters (pij, the nugget parameter Q, and the ith element 
of 0, 9i. For each jth principal component, we estimate the parameters Kj, <p\j ,..., <p ( j 0 , and (j 
by maximizing the log likelihood function t(wj(Q\), ..., Wj(O p )\Kj, 4>\j, ..., 4>dj,(j) given by the 


Gaussian process for Wj(0\),... ,Wj(O p ) defined by the covariance function in (14). Here we use 
the exponential covariance function because we expect that there is not enough information in 
the binary data to infer the smoothness parameter in Matern. Our choice of smoothness pa¬ 
rameter 0.5 for the Matern covariance function (therefore an exponential covariance function) is 
a conservative one. We have also conducted a leave-10-percent-out cross validation experiment 
and found that an emulator with the exponential covariance function shows a better interpola¬ 
tion performance than the squared exponential covariance function. The computational cost for 
each likelihood evaluation is flops, which is 4 x 499 3 « 4.14 x 10' flops in our application. 

The Jy-dimensional emulator process r/ (0 , W) for the logistic PC scores at any new param¬ 
eter value 0 is obtained as usual from the conditional distribution. Then the n-dimensional 
emulator process for the natural parameters at 0 is 


1 (0) = n + K y r 1 (0,W), 


(15) 


and the corresponding probability for the predicted model output at each location to be one is 


Pj (0) = P(Y(0, Sj ) = 1) = 


exp ( 7 j(0)) 

1 + exp (jj(0)) ’ 


where 7 j(0) is the jth element of "y(0) in (15) 
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4.3 Calibration Using Logistic PCs 

Once the emulator for the logistic PCs is available, we set up the following calibration model for 
the n-dimensional vector of the natural parameters A for observational data using the emulator 

H + Kyri(0,W): 

X = ti + K y r] + d, (16) 

where 6 is the best fit value for the natural parameters for observational data, r/ = rj (6, W) is 
the emulated natural parameter vectors at the value of input parameter 6, S is the discrepancy 
term that represents structural error between the natural parameters for the model output and 
those for the observational data. The intercept /i and the basis matrix K. y are given from the 
emulation stage. To get around the challenges in integrating out <5 described above, we use a 
basis representation for the discrepancy term such that 


<5 = K d v, 

with the nxJ d basis matrix and the J^-dimensional random coefficient vector v r\_/ N(0,ajl Jd ) 
Similar to the independence assumption between the emulated vector and the discrepancy term 
made for the model in Q, we assume that (v,<r^) is independent of ( rj,G,W ). 

The posterior density given the observational data Z (see Supplementary Section A1 for 
details about its derivation) can be written as 


vr(?7, 6. v, cxrf|Z, W) oc L (Z\rj, v) / (rj\0, W) / (v|<r d ) f{0)f(p d ), 


(17) 


where the likelihood function L (Z| rj,v) is 


L(Z|t7,v) oc J j 


exp(Aj) \ z(Sj) f 1 V _Z(Bi) 


L V1 + exp(Aj) J V1 + e x p(Ay) 


with A j being the jth element of A = /x + K y r) + 5 , and / {rj\6, W) and / (v|crd) are 


; (r ' 19 ' w) K j} 77w exp 

Jd ^ 

f iy\(T d ) ocTT —exp 

- LJ - (Td 

1=1 a 


2 °dJ ' 
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where n ri .k(0) and crj fc (0) are the mean and the variance of the kill element of ri(0, W), vi is 
the Zth element of v, and crj is the common variance for v\.. . vj d . For the priors we assume 
that f(9) is a (/-dimensional uniform density whose support covers all plausible values of 6 
and f((Td) is an inverse-Gamma density with the shape parameter ad and the scale parameter 
bd . The functions L(Z\rj,v), f (rj\9,W), and f (v\ad) are based on proper densities, hence 
the posterior is proper as long as the priors f{9 ) and f(<Jd) are proper. The observations 
Z(s\Z(s n ) are conditionally independent given the values of the random coefficients rj 
and v and the dependence between them is accounted for through the basis matrices K y and 
K d- Based on the specifed posterior one can infer the parameters via MCMC. 

Finding the basis that can precisely represent the discrepancy process with a small 
number of basis vectors is important for computational tractability. To this end, we use the 
following procedure that is inspired by the definition of the discrepancy term: 


1. For each location s j, compute the signed proportion of mismatch between the model out¬ 
put and the observational data as rj = Y^=i sgn(y(0j, s j) — Z{sj))I(Y(6i , s j) / Z(sj))/p 
where sgn is the sign function. 

2 . Set Jd = 1 and hence is reduced to a re-dimensional vector k^. Define the jth element 
of as log(j^^) if | rj | > c, or as 0 if \rj\ < c. The random coefficient vector v is also 
reduced to a univariate random variable v ~ N(0,a^). 


The goal of step 1 is to capture a common discrepancy pattern that is persistent across all 
parameter settings. For this procedure to work the design points 0i,...,6 p have to be a 
representative sample of plausible values of 6, and we expect that this holds for the Latin 
hypercube design that we are using here. In step 2 we apply logit transformation log f° r 

—1 < rj < 1 to translate the identified discrepancy pattern into logit scale since the discrepancy 
term k dV needs to be specified for the natural parameters. Note that the cut-off value c in step 
2 needs to be determined a priori. Choosing a value of c that is too large results in loss of 
important discrepancy patterns in k^ and, on the other hand, a value of c that is too small 
causes identifiability issues in estimating 6. In the simulated examples described in Section 


5.1, we have found that the value of c = 0.5 strikes a nice balance, and yields a discrepancy 


basis that capture the true simulated discrepancy patterns reasonably well (See Figure A1 in 
the supplement). This value of c is a heuristic but it appears to work well for the variety of 
examples to which we have applied our methods. We set c = 0.5 for the remainder of this 
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paper. Note also that log( \ ) can take a positive or negative infinite when rj = 1 or rj = — 1. 

To prevent this, we set rj = 0.999 x sgn(rj), where sgn(rj) is the sign of rj , when computing 
l°g(l3fr) ^ \ r j\ is greater than 0.999. 


4.4 Procedure Summary 

Let Y be an n x p matrix of the computer model output where the columns represent dif¬ 
ferent spatial locations and the rows represent different parameter settings and Z denotes an 
n-dimensional vector of observational data (defined in Section [3]). We summarize the overall 
emulation-calibration procedure described in this section as follows: 


1. Dimension-reduction of computer model output: Compute logistic principal components 


of Y to reduce the dimensionality from n to J y (details in Section 4.1) by maximizing the 
likelihood function in ([7]). We store the principal component scores in apx J y matrix 
W, the principal component basis in an n x J y matrix K^, and the location means in an 
n-dimensional vector /i (see ([6]) for their definition). 


2. Gaussian process-based univariate emulators: We construct a one-dimensional Gaussian 


process emulator for each column of W (details in Section 4.2) with the covariance function 


in (14). Denote the ,7,,-dimensional vector of the emulated principal components at a new 


parameter value 6 by r](6 , W). 


3. Constructing a discrepancy basis: Find the discrepancy basis vector by comparing each 


row of Y and Z (details in Section 4.3). Denote the coefficient for the discrepancy basis 
vector by v. 


4. Calibration using Markov chain Monte Carlo: Based on the posterior density in ( |17[ ) infer 
the input parameters 0 , the principal components rj, the discrepancy coefficient v, and 
the discrepancy variance a y using the Metropolis-Hastings algorithm (details in Section 

O). 


Our procedure requires estimating only J y + d + 2 parameters in the calibration step. For our 
ice sheet model calibration problem, this corresponds to only 10 + 10 + 2 = 22 parameters. 
Note that it would require estimating n + d + 1 = 3,182 + 10 + 1 = 3,193 parameters without 
dimension reduction. 
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5 Implementation Details and Results 


In this section we study the application of our approach to simulated examples and then calibrate 
the PSU3D-ICE model based on the real observational data. We provide descriptions for each 
of the calibrated parameters as well as scientific interpretations of our results. 

5.1 Simulated Examples 

We now provide a discussion of our investigation of our approach in the context of simulated 
examples. We describe our construction of simulated examples and results from the applications 
of our methods to these examples, along with implementation details for our method. We apply 
our method to (i) a simple synthetic ice sheet model example to study the properties of our 
method, and to (ii) the PSU3D-ICE to estimate the input parameters of the model and make 
projections for future ice sheet volume changes in the Amundsen basin region. 

Simple Numerical Example. We begin with a simple numerical example. Unlike the real ice 
sheet model calibration problem, this example does not have complicating issues such as lack 
of identifiability of input parameters or sparsity of design points in the parameter space, and 
hence it provides a nice first test for our calibration methods. 

We defined the model output as Y (0, s) = I ( U(9 , s)) where 

U(0,s) = \J 1 — s'l — — Oi + 6 * 2(52 + 1.5), 

with s = (si, S 2 ) and 6 = ( 6 \, 62 ) for 0.3 < 6 \ < 0.65, 0 < 62 < 0.2, 0 < si < 1, and —1.5 < S 2 < 
1.5. We form an ensemble of model outputs with 100 runs where the design points for (6fi, O 2 ) 
are on a regular 10 x 10 lattice covering the entire parameter space 0 = [0.3,0.65] x [0,0.2], 
Each model run Y ( 9,) is a spatial pattern on a 30 x 30 regular lattice that covers the entire 
spatial field S = [—1,1] x [—1.5,1.5] and hence n = 900. See Figure A2 in the Supplement for 
16 sample model runs. The first parameter 6 \ controls the overall size of the grey area and the 
second parameter 62 changes both the shape and the size of the grey area. We set the number 
of principal components to J y = 10 based on leave-10-percent-out cross-validation; using more 
than 10 principal components does not significantly improve the emulation performance in our 
cross-validation experiments. The overall misclassification rate is less than 5%. We construct 
a synthetic observational data set by choosing the model run with (Q \, 62 ) = (0.494,0.089) as 
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the synthetic truth and contaminating the corresponding model output with noise generated 
from the discrepancy process, defined as a zero-mean Gaussian process with an exponential 
covariance function with the partial sill of 1.5, the range of 0.03 and the nugget of 0.00001. 
This affects the resulting binary pattern considerably by changing roughly 10% of the original 
output. We repeat the same experiment with 10 different realizations from the discrepancy 
process and the results are summarized in Figure A3 in the Supplement. The results show 
that, on average, we can recover the true input parameter setting reasonably well in this simple 
numerical setting. 

Realistic Example Using the Ice Sheet Model PSU3D-ICE. We now work with a more realistic 
synthetic observational data set to verify that our calibration approach works well in the context 
of the real ice sheet model calibration problem. Here we use the model runs from the PSU3D- 
ICE model described in Section [2] to generate the example. This estimation problem is more 
challenging than the example described in the previous section because it involves complicated 
interactions between the input parameters. Also, the input parameter settings are much more 
sparsely distributed in the parameter space (499 ensemble members in 10 dimensional space), 
which potentially makes emulation, as well as calibration, more difficult. 

We select model run #394 as the synthetic truth, because it is among the model runs closest 
to the observational data. We count the number of pixels where the observational data and 
the model output have the same value, i.e. we compute = T(0j,s j)) for each 6 , 

where /(.) is the indicator function. A model run is close to the observations if this count is 
large. Since in reality the model-observation discrepancy operates in terms of the actual ice 
thickness patterns (the binary ice-no ice patterns are derived from these data), we construct 
the discrepancy pattern in terms of ice thickness and translate it into a binary spatial pattern. 
We use the following steps to construct a realistic model-observation discrepancy. 

1. By comparing the original model runs and the observational data that are spatial patterns 
of ice thickness, find the top 60% model runs with the smallest root mean squared error 
(RMSE) values. 

2. For each spatial location, take the average of the difference between the observational 
data and the model runs selected in step 1 to find the common discrepancy pattern in ice 
thickness. 
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3. Subtract common discrepancy pattern found in step 2 from the ice thickness output cor¬ 
responding to the synthetic truth. The resulting pattern is synthetic observational data 
for ice thickness. 

4. Dichotomize the ice thickness pattern computed in step 3 at 0 to obtain synthetic obser¬ 
vational data for ice-no ice pattern. 


The procedure described in Section [4.3| recovers the superimposed discrepancy reasonably well 
(Figure A1 in the Supplement, lower panels). 

We choose to use 10 principal components for dimension reduction based on the same cross- 
validation criteria described in the previous section. The leave-10-percent-out cross-validation 
confirms that our emulator can predict the model output precisely (Figure [2]), with an overall 
error rate of 4%. To translate the probabilities pi(6),... ,p n (6) from our emulator into predic¬ 
tions for the binary outcome, we simply dichotomize them using 0.5 as a threshold. We give a 
non-informative prior IG( 2,3) to crj and flat priors to the input parameters in 0. The range 
of each flat prior is restricted to the range for the design points 0 1 ,..., 0 p . Using the standard 
Metropolis-Hastings algorithm we generate MCMC sample for 6 while integrating out rp v, and 
(jj. For each experiment we obtain MCMC chain with 1,000,000 iterations. Looking at MCMC 


standard errors (Jones et al. 2006 Flegal et al. 2008) suggests that the sample size was large 
enough to estimate posterior means while accounting for autocorrelations. We also overlaid the 
marginal density plots obtained after half the run over the marginal density plots obtained after 
the entire run, confirming that the plots do not show significant changes and the results are 
therefore stable. The computation takes about 120 hours for the entire MCMC run on a system 
with Intel Xeon E5450 Quad-Core 3.0 GHz.. 

The pairwise joint densities of input parameters are shown in Figure [3j Since the scientific 
interest is mainly in the four parameters, ocean sub-ice-shelf melt factor (OCFAC), calving 
factor (CALV), basal sliding coefficient (CRH), and asthenospheric relaxation e-folding time 


(TAU), we display the results for only these parameters (see Section 5.2 below for the reason). 
The density plots indicate that the marginal posterior densities for the first three densities peak 
around the true parameter values, but the marginal densities for TAU do not recover the true 
values and show large dispersion comparing to the other parameters. We hypothesize that the 
binary spatial patterns are not informative about TAU. The probable reason for parameter TAU 
not being well constrained is that the statistical analysis above is limited to calibration with the 
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modern observed ice state. TAU is the e-folding time scale of asthenospheric relaxation below 
the varying ice load (see Table [I] and Section 5.2 below), and its primary effect on ice evolution 
occurs during the main deglacial retreat phase, 12 to 8 ka. The ice retreat is nearly complete 
by 8 ka in most of our model runs, leaving ample time for the asthenosphere to relax towards 
the modern state, regardless of the exact value of TAU. Calibration and statistical techniques 
using past grounding-line positions vs. those estimated from geologic data over the last 15 kyr 
(Larter et al. 2014) will likely be able to diagnose and constrain this important parameter. 

Using the MCMC chain for the input parameters, we generate projections for WAIS volume 
change. Here we present projection results for 500 and 1000 years from the present day. We 
build an emulator for ice volume changes generated from the PSU3D-ICE and convert the 
MCMC chain for input parameters into projections using the emulator. The results in Figure 
[4] show the predictive distribution of the ice volume change projections corresponding to the 
estimated posterior density. For comparison, we have also computed the predictive density by 
assuming that all the parameter settings in the model ensemble are equally likely. We compare 
these densities with the assumed true ice volume change projection values produced by the 
parameter settings for the synthetic truth. For both time points (500 and 1000 years from the 
present) the 95% prediction intervals cover the synthetic true ice volume change projections. 
The mode of the predictive density matches well with the true ice volume change projection 
for 500 years, but the mode for the ice volume change for 1000 years does not match with the 
assumed true value. We hypothesize that this is due to the large uncertainty for TAU, which 
controls the long term ice sheet behavior. We provide a longer discussion of this in the next 
section. 


5.2 Calibration of PSU3D-ICE Using Observational Data 


We now calibrate PSU3D-ICE using the observational data. The calibration results based on 
the observations from the Bedmap2 dataset (Fretwell et ah, 2013) are shown in Figure [6] and [Tj 


Parameter Description. Out of the 10 ice-sheet model parameters that were varied and 
inferred we focus on four important parameters. These parameters are highly uncertain even 
though they have strong effects on model behavior. Most parameters and inputs that concern 
interior grounded-ice processes, inland of the modern grounding zone, are relatively well known, 
based on modern data and well established in previous model-data studies. Examples are 
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ice rheology, surface accumulation, and basal sliding through inverse procedures. However, 
processes that affect the floating ice shelves, and also properties of past grounded ice that 
advanced out over the Antarctic continental shelf at glacial maxima, are much less well known, 
and we focus on important parameters in this category. The 4 parameters and processes are: 


Ocean sub-ice-shelf melt factor (OCFAC): Oceanic melting at the base of floating ice 
shelves, due to warm waters flowing from the open ocean into the cavities below the ice 
shelves. This affects ice shelf thickness, and changes buttressing of interior ice. Recent 
increases in ocean melt rates are considered to be the main cause of drastic ongoing ice 


retreat in the ASE sector of WAIS (Pritchard et al., 2012; Rignot and Jacobs, 2002). For 
small (large) OCFAC values, oceanic melting is reduced (increased), ice shelves thicken 
(thin), discharge of interior ice across the grounding line decreases (increases), and ground¬ 
ing lines tend to advance (retreat). 


• Calving factor (CALV): Calving of icebergs at the oceanic edge of floating shelves. This 
process is not well understood, probably depending on pre-existing fracture regimes, large- 
scale stresses, and hydrologic conditions, yet has important effects on ice-shelf extent and 
feedback on buttressing and internal ice, both for the modern state and also for past 
and future variations. There is little consensus on calving parameterizations. We use 
a common approach based on (highly parameterized) crevasse depths and their ratio to 
ice thickness (Benn et al. 2007; Nick et ah, 2010). For small (large) CALV, calving is 
decreased (increased), producing more (less) extensive floating shelves, and greater (lesser) 
buttressing of interior ice. 


Basal sliding coefficient (CRH): Basal sliding of grounded ice at the interface with bedrock. 
For beds at the ice pressure-melt point, for a given driving stress, sliding depends mainly on 
whether the bed is hard and non-deformable (crystalline rock, sticky, slow sliding) or soft 
and deformable (unconsolidated sediment or till, slippery, fast sliding). The coefficient 
values in the basal sliding law relating basal stress to velocity have strong effects on 
the overall ice profile shape and thicknesses, and affect the overall sensitivity to climate 


change. Coefficients under modern grounded ice are deduced by inverse methods (Pollard 


and DeConto, 2012a Morlighem et ah, 2013), but they are relatively unconstrained for 


modern oceanic beds, where grounded ice advanced during the last glacial maximum 20 
to 15 ka. Most oceanic areas around Antarctica are covered in deformable sediment today, 
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due to marine sedimentation and to generation and transport by previous ice advances. 
For these regions, there is a wide plausible range of relatively large sliding coefficients that 
have significant effects on model results. Here, we vary the sliding coefficient CRH just 
for modern oceanic areas. 


Asthenospheric relaxation e-folding time (TAU): The ice sheet evolution on 10 3 year 
timescales depends strongly on the bedrock response to varying ice load. During deglacial 
retreat, the bedrock rebounds upwards due to the retreating ice load, but with an O(10 3 ) 
year lag, and this lag has an important positive feedback via deeper bedrock and grounding¬ 
line depths. As in many ice models, this bedrock response is represented by a simple Earth 
model consisting of an elastic plate (lithosphere) over a local e-folding lagged relaxation 
towards isostatic equilibrium (asthenosphere). Based on more sophisticated Earth mod¬ 


els, the asthenospheric e-folding time scale is commonly set to 3 kyr (e.g., Gomez et al. 


2013), but note that recent studies suggest considerably shorter time scales for some West 


Antarctic regions (Whitehouse et al. 2012a). 


Results. The estimated discrepancy pattern by the procedure described in Section |4.3| is 
shown in Figure [5j Much of the discrepancy pattern in Figure [5] is due to a known shortcoming 
of the ice model that tends to fill narrow bays and passages between mainland and nearby 
islands with grounded ice, where there is only floating ice in the real world. This is partly due 
to the coarse model grid size (20 km) not sufficiently resolving the passages, and also presumably 
due to the lack of fine-scale variations in oceanic basal melting that could well exist in these 
passages due to local ocean currents. 

As seen in the pairwise density plots in Figure [6j the best values of ocean melt coefficient 
(OCFAC) and calving factor (CALY) are roughly in their mid-range, i.e., quite close to the 


values in the nominal model that relied on prior tuning by the model developers (Pollard and 


DeConto, 2012a). There is a slight indication that OCFAC should be somewhat smaller than 


the value found by Pollard and DeConto (2012a). The best CALV value, 0.5, corresponds to 
an unchanged nominal value. 

The bedrock (asthenospheric) e-folding rebound time TAU is not constrained well in Figure 
[6j Although the uppermost values (corresponding to 6 to 7 kyr in the original scale) are 
excluded, values below 6 kyr have similar probabilities. As discussed above, this is expected 
because TAU primarily affects past ice distribution during deglaciation, with little effect on 
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modern ice. There is a slight tendency to favor small TAU (less than 0.2, corresponding to less 


than 2 kyr), compared to the nominal prior value of 3 kyr (Pollard and DeConto, 2012a), but 
this is a weak tendency. 

There is one strong signal in Figure[6| the best-fit basal sliding coefficient CRH for continental 
shelf (modern ocean) regions, is found to be quite high, 0.8, corresponding to 10~ 6 nr a -1 P a -2 . 
The narrowing of the uncertainty for CRH is valuable information for models. Even though the 
validation here is based only on results for the modern ice sheet (for which there is no direct 
effect of CRH if the results are realistic, because all modern ice should be floating in the affected 
regions), CRH still has a significant influence on the modern state by affecting ice profiles during 
the recent evolution over the last few thousand years. The relatively high (slippery) value of 
10 -6 is consistent with geologic data at LGM, indicating that grounded ice thicknesses across 
the modern Ross Sea Embayment were relatively thin at LGM (deduced from geologic data 
in the Transantarctic Mountains and other outcrops; Ackert Jr. et al.[ |2013| Anderson et ah 


2013). 


In the ice-volume change projections for 500 and 1000 years from present (Figure[7]), the left- 
hand tail of lower ice volumes may seem counter-intuitive, given the imposed climatic warming 
after year 0 (present day). This tail is probably due to a subset of highly unrealistic runs, 
whose grounding lines retreat drastically inland beyond the modern coastline much too early, 
100 ’s to 1000’s years before present, and then re-advance due to the lagged bedrock rebound 
(discussed above), producing modern grounding lines close to observed despite very unrealistic 
past behavior. This likely happens in runs with the longest TAU time-scale values, 7 kyrs. As 
mentioned above, calibration based on past geologic data would reduce this tail. 

The difference between the mode at 1000 years from present and the assumed-true value 
(Figure [4]) could well be due to the unconstrained TAU values in Figure [3] and [gJ As well as 
affecting past behavior, TAU also has an effect on long-term future retreat via the feedback 
on grounding-line depths discussed above. Nevertheless, the basic projection of drastic future 
retreat and a few meters (less than 10 meters) equivalent sea-level rise within the next 1000 
years is a robust result, and occurs in nearly all runs with reasonably realistic past behavior. 
This is consistent with other recent studies of future ice-sheet retreat in the Amundsen Sea 


Embayment sector (Joughin et ah, 2014 Rignot et ah, 2014) 


Caveats. Similar to other calibration frameworks, care is needed in specifying the model- 
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observation discrepancy to avoid identifiability issues. In our calibration approach specifying 
the discrepancy model requires a pre-determined tuning parameter c, which determines whether 


a specific pixel exhibits a persistent model-observation discrepancy (see Section 4.3). While our 
choice of c = 0.5 appears to be somewhat ad-hoc, we have verified through a detailed study of 
multiple simulated examples that this choice of c generally yields a discrepancy pattern that 
captures well the true discrepancy pattern. In addition to the simulated example (Section |5.1[ ) 
and the PSU3D-ICE model used in this paper, we have applied our approach with c set to 0.5 


in the context of a completely different ice sheet model (SICOPOLIS; Greve, 1997; Greve et al. 


2011 ) and verified that the identified discrepancy is similar to the true discrepancy pattern used 


to create the simulated examples. 

Our scientific conclusions are subject to the following usual caveats. It should be emphasized 
that the extensions of our runs beyond present are projections based only on a single, simpli¬ 
fied greenhouse-gas scenario. Various greenhouse-gas scenarios for the future are available as 


Representative Concentration Pathways (RCPs Moss et al. 2010 Meinshausen et al. 2011), 
which can be used with climate models to provide forcing for the ice sheet model. However, the 
RCPs themselves range over a wide envelope (RCP 2.6, RCP 4.5, RCP 6.0, RCP 8.5) depending 
on future societal activity, from strong through intermediate mitigation to “business-as-usual”. 
These scenarios should ideally also be treated as sources of additional projection uncertainty. 
This is a complex additional challenge and a potentially important direction for future work. 
In this preliminary study we limit the variable parameters to physical model attributes, and 
specify a relatively simple future climate warming. Note also that the volumes in Figures [4] and 
[T] are only for the WAIS domain of these nested model runs, not for all Antarctica (although 
the East Antarctic contribution is expected to be comparatively small; Pollard and DeConto, 


2009). 


Our calibration results are based only on the modern day grounding-line geometry and 
it might be possible to improve the results by using additional information such as the past 
grounding line positions. We expect that this will lead to a better constrained probability den¬ 
sity of TAU. As we mentioned above, one future direction is to formulate a rigorous calibration 
approach that enables us to combine these multiple sources of information. In addition, we 
use only the grounding-line geometry and hence ignore the spatial patterns of floating ice (ice 
shelves) and the distribution of ice thickness. Incorporating these additional sources of infor¬ 
mation requires calibration using other types of non-Gaussian data such as bounded spatial 
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data with many zeros and spatial data with multinomial responses; this is beyond the scope 
of this paper. Also, methods combining calibration with other types of past geologic data will 
be challenging but necessary for comprehensive model validation (mainly Relative Sea Level 
records, GPS uplift rates, and cosmogenic age-elevation data;|Briggs and Tarasov, 2013, Briggs| 
et ah, 2013, 2014 Whitehouse et ah, 2012b|a ). 


6 Summary 

We have formulated a novel computer model emulation and calibration approach for high¬ 
dimensional binary spatial data. We have applied our approach to calibrate a computer model 
that describes the dynamics of the West Antarctic ice sheet. The results of this calibration 
provide some understanding of the role of important parameters in the model, and also allow 
us to produce preliminary projections, with uncertainties, of future ice sheet volume change. 
Using a logistic principal component analysis and basis representation, we can efficiently han¬ 
dle high-dimensional binary spatial data; our dimension-reduced approach therefore addresses 
challenges in handling high-dimensional binary spatial data by removing the need to simulate 
underlying high-dimensional processes and avoiding computations with large matrices. The lo¬ 
gistic principal component analysis is a special case of principal component analysis for the one 
parameter exponential family and can therefore be easily extended to other types of spatial data 
such as spatial count processes. Incorporating complicated model-observation discrepancies is a 
challenging problem in general. We allow for a flexible but computationally expedient model for 
discrepancy. Our study of simulated examples shows that our approach can estimate the input 
parameters reasonably well even in the presence of considerable, realistic model-observation 
discrepancies. The application of our methods to real ice sheet observational data allows us to 
characterize uncertainties in both our model parameter estimates as well as our projections for 
future ice volume change. 
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Parameter 

Definition 

OCFAC 

CALV 

Ocean sub-ice-shelf melt factor (.1 to 10) (non-dim). K in Eq. 17, PD12. 
Calving factor (.316 to 3.16) (non-dim), multiplies combined crevasse depth to 

CRH 

thickness ratio r in Eq. B7, PD15. 

Basal sliding coeff. modern ocean bed (10~ 9 to 10~' 5 ) (m yr -1 P a -2 ) C in Eq. 

11, PD12. 

CRHASE 

Extra multiplicative factorxCRH, inshore ASE (10~ 3 to 1) (non dim), multiplies 
C, only in inshore Amundsen Sea Embayment, Eq. 11, PD 12. 

TAU 

LITH 

GEO 

SUBPIN 

Asthenospheric relaxation e-folding time scale= (1 to 6 kyr). t in Eq. 33, PD12. 
Lithospheric stiffness (10 23 to 10 25 ) (N m). D in Eq. 30, PD12. 

Geothermal heat flux (50 to 90) (mW m~ 2 ). see section 2.7, PD12. 

Sub-grid ice-shelf pinning factor (0 to 4) (non-dim), multiplies Sdev hi Eq. 13, 
PD12. 

USCH 

LAPSE 

Grounding-line flux factor (0.5 to 5) (non-dim), multiplies q q in Eq. 8, PD12. 

7 (lower-case greek) in Eq. 34a, PD12. atmospheric lapse rate (-.005 to -.010) 
(°C m” 1 ). 


Table 1: Ice-sheet model parameters varied in the LHC ensemble described in Section [ 2 } The 
parameter ranges and units are shown in parentheses. The parameter values are uniformly 
spaced within each range, but for some variables, they are log-linear (i.e., the log of the param¬ 
eter is varied uniformly), for OCFAC, CALV, CRH, CRHASE, LITH, and USCH. Parameters 
called ’’factors” are non-dimensional, multiplying existing terms. The corresponding symbols 
and equation numbers in Pollard and DeConto (2012a) and Pollard et al. (2015) called PD12 
and PD15 respectively, are also shown in the table. 
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Observational Data Model Output from Run No.394 



Figure 1: Spatial patterns of grounded ice from the observational dataset 
(left) and an example PSU3D-ICE run (right). The gray pixels represent 
The model appears to be capable of reproducing observational patterns. 


(Fretwell et al., 2013) 


ice-covered locations. 
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Model Output from Run No.67 


Emulated Output for Run No.67 



Model Output from Run No.491 



Emulated Output for Run No.491 



Figure 2: Two examples of leave-10-percent-out cross-validation results to study the perfor¬ 
mance of our logistic PCA based emulator. The gray pixels show ice-covered locations, (a) 
Comparison between the original (left) and the emulated (right) ice coverage patterns for model 
run no. 67. (b) The same comparison for model run no. 491. The graphical comparison shows 
that the emulated model output (right panel) approximates the original model output (left 
panel) fairly closely. Similar results hold for the other model runs. 
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2-D Posterior Densities for Input Parameters 
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Figure 3: The pairwise joint posterior densities for the input parameters from the calibration 
results for the simulated example in Section 5.1 We plot only the following four parameters that 


we are most interested in: ocean sub-ice-shelf melt factor (OCFAC), calving factor (CALV), 
basal sliding coefficient (CRH), and asthenospheric relaxation e-folding time (TAU) (see Section 
5.2 for parameter descriptions). The white dashed lines show the input parameter settings 
for the synthetic truth. The marginal densities peak around the assumed true values for the 
parameters. The lone exception is TAU, indicating that the modern ice coverage pattern is not 
informative about this parameter. 
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Projections based on synthetic data (500 years) 
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Projections based on synthetic data (1000 years) 



Figure 4: The ice volume change projections for 500 years (left panel) and 1000 years (right 
panel) from the simulated example described in Section |5T| based on the calibration results (solid 
red lines) in Figure [3] and no calibration (dashed green lines). The vertical black line shows the 
true projection value generated by the synthetic truth. The parentheses on the horizontal bars 
represent 95% prediction intervals and the vertical bars the ranges. The slight bias in the mode 
of the predictive distribution for the volume changes for 1000 years is likely due to the bias for 
the posterior mode of TAU (see Figure [3]). 
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Figure 5: Estimated discrepancy between PSU3D-ICE model and the Bedmap2 observational 


dataset (Fretwell et al., 2013) by the procedure described in Section 4.3 
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2-D Posterior Densities for Input Parameters 
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Figure 6: The pairwise joint posterior densities for the input parameters from the calibration 
results based on the real observational data described in Section [2j Again we plot only the 
following four important parameters: ocean sub-ice-shelf melt factor (OCFAC), calving factor 
(CALV), basal sliding coefficient (CRH), and asthenospheric relaxation e-folding time (TAU) 
(see Section 5.2 for parameter descriptions). Similar to the results in Figure [3j the posterior 
densities have tightly constrained high density regions except for TAU. 
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Projections based on observations (500 years) 



Projections based on observations (1000 years) 
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Figure 7: The ice volume change projections for 500 years (left panel) and 1000 years (right 
panel) based on the calibration results (solid red lines) using the real observational data shown 
in Figure [o] and no calibration (dashed green lines). The parentheses on the horizontal bars 
represent 95% prediction intervals and the vertical bars represent the ranges. 
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